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FULL WAVE CALCULATIONS OF 
THERMOSPHERIC NEUTRAL AIR MOTIONS 

H. Volland 
ABSTRACT 

Full wave calculations have been performed within the frequency range of 
gravity waves (lO -4 < co < 10~ 2 sec -2 ) for a thermospheric model between 120 
and 700 km altitude. In this altitude range gravity waves are coupled with heat 
conduction waves. Reflection, transmission, conversion and coupling from one 
wave type into the other one is described by the elements of the scattering 
matrix. The dependence of these elements on frequency and angle of incidence 
is discussed in full detail. The transmission coefficients of gravity waves cal- 
culated by full wave theory are compared with simple ray calculations and show 
that ray treatment is an excellent approximation for obliquely upward propa- 
gating gravity waves and that gravity waves predominate throughout the thermo- 
sphere. The thermosphere reacts like a selective filter with respect to upward 
propagating gravity waves with optimal transmission at k x ~ 2to (a> = angular 
frequency; k x = horizontal wave number). Vertically downward propagating 
heat conduction waves are mainly converted into upward reflected gravity waves 
(at higher frequencies) or mainly coupled into downward propagating gravity 
waves (at lower frequencies). 
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FULL WAVE CALCULATIONS OF 
THERMOSPHERIC NEUTRAL AIR MOTIONS 

1. Introduction 

Satellite drag measurements of the neutral air density (Jacchia, 1959; 
Priester, 1959) have revealed for the first time that large scale diurnal varia- 
tions of the density in thermospheric heights between 150 and at least 1200 km 
altitude exist, thermally driven by EUV heating from the sun (Harris and 
Priester, 1962). Hines (1960) put forward the hypothesis that free internal 
gravity waves of smaller periods generated within the lower atmosphere can 
propagate into the thermosphere. He suggested that travelling ionospheric dis- 
turbances observed within the F -region (see e.g., Heisler, 1967) are the response 
of the ionospheric plasma to these neutral air waves. Theoretical calculations 
especially those dealing with the ionospheric effect of the Russian H-bomb ex- 
plosion of October 30, 1961 (Kohl, 1964; Row, 1967; Hines, 1967a) show sur- 
prising good agreement between the observed variations of the F -layer critical 
frequency at different stations and the theoretically expected dispersion effect 
of internal gravity waves depending on angle of incidence. Hines (1967a) could 
even explain the apparent period of the disturbance by a simple ray approach of 
the direction of energy propagation of gravity waves. 

Recently, radar backscatter measurements (Thome, 1964) and HF-Doppler- 
effect measurements (Georges, 1967) as well as neutral air density observations 
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by pressure gauge measurements on board of Explorer 32 (Newton et al. 1968) 
have shown that wave like structure in the thermosphere reaches altitudes as 
high as 600 km. Hines (1967b) and Newton et al. (1968) using a ray treatment ex- 
plained these waves again as free internal gravity waves propagating from the 
lower atmosphere up into the thermosphere. The most striking feature of these 
large altitude waves is the bending of the wave normals toward a horizontal 
direction with increasing altitude. In a ray treatment this bending is a dispersion 
effect due to energy dissipation of the waves [Hines (1967b); Newton et al. (1968)] . 

From the work of Harris and Priester (1962) and of Pitteway and Hines (1963) 
it can be suggested that heat conductivity is of predominating influence for the 
energy dissipation of neutral air waves. It can be shown numerically that the 
coefficient of viscosity can be neglected as compared with heat conductivity in 
thermospheric heights below 600 km and within the frequency range of gravity 
waves (Volland, 1968). 

Under the influence of heat conduction, two different pairs of wave modes 
exist within the thermosphere: acoustic -gravity waves and heat conduction 
waves. They are coupled with each other at any height because their eigenvalues 
depend on the ratio 

K 

P 

(k = coefficient of heat conduction; P = mean pressure) 
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which increases with height, and therefore the thermosphere behaves like an 
inhomogeneous medium. Heat conduction waves are evanescent waves below 
about 200 km altitude. At higher altitudes, however, their attenuation rate is of 
the same order of magnitude as the attenuation rate of gravity waves. For these 
two reasons, the response of the thermosphere to the different wave modes in- 
cluding coupling and reflection processes must be determined by a full wave 
treatment, and a ray approximation has to be justified by these exact calculations. 

Full wave calculations of neutral atmospheric waves in altitudes below 200 km 
have been performed by Midgley and Liehmon (1966). They calculated reflection 
coefficients of gravity waves and eliminated the evanescent heat conduction waves 
and viscosity waves by a special mathematical technique. 

As pointed out above, this technique, which is fully justified at altitudes 
below 200 km, is not acceptable at greater heights because there heat conduction 
waves are no longer of an evanescent type. 

This paper deals with a full wave treatment of neutral atmospheric waves at 
thermospheric heights taking into account heat conduction but neglecting viscosity, 
ion drag and Coriolis force. The negligence of the coefficient of viscosity is 
allowed at heights below 600 km and within the frequency range of gravity waves 
(o> £ IC 2 sec ' 1 ). The negligence of ion drag and Coriolis force is possible at 
frequencies 
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where ft is the rotation period of the earth and v is the collision frequency be- 
tween neutrals and ions. 

Our numerical calculations are extended to an altitude range of 120 I Z | 
700 km and a frequency range of 10 4 < « < 10 2 sec 1 . We shall compare 
these full wave calculations with simple ray approximation and shall mark the 
ranges of frequency, angle of incidence and wave type where a ray treatment is a 
sufficient approximation for neutral air wave propagation. 

The mathematical background of full wave calculations is presented in 
Appendices A to E. For an understanding of the numerical results given in 
Section 4 only the contents of Sections 2 and 3, are necessary. 

2. The Elements of the Scattering Matrix 


As shown in Appendix A to C two pairs of characteristic waves exist within 
the neutral thermosphere which we collect by a column matrix [see Equation (C-l)] . 


/ a G \ 
a H ' 
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where a Q , b G are ascending and descending gravity waves respectively, and a H , 
b H are ascending and descending heat conduction waves , respectively. 

Our model atmosphere consists of three regions (see Figure 1). The regions 
below Zj (region I) and above Z n (region III) are considered as homogeneous 
(G = const) and isothermal (T = const) atmospheres extended infinitely into the 
half spaces. The region between Z I and Z u (region II) approximates part of the 
real atmosphere. Its parameters approach the parameters of regions I and III 
at the boundaries Zj and Z n , respectively. 


Coupling between the four characteristic waves then only occurs within 
region II, while in regions I and in the waves are uncoupled with each other. 

The connection between the four waves at the lower boundary Z x with the four 
waves at the upper boundary Z n of region n is given by Equation (D-6): 

a 11 (Z u ) = Pj 11 a 1 (Zj) (2) 


We are mainly interested in reflection and transmission coefficients of the 
different waves. These coefficients are the elements of the scattering matrix 
well known in electromagnetic wave propagation (see e.g., Volland, 1962). The 
scattering matrix connects the waves going into region II with the waves coming 
out of region It: 


b 1 = R 1 o 1 + T n b 11 
0 n = t 1 a 1 + R 11 b 11 


(3) 
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where we introduced column matrices of ascending and descending waves: 




R 1 and R 11 are (2 x 2) — reflection matrices, and T 1 and T 11 are (2x2) — 
transmission matrices . We collect these last four matrices into a (4 x 4) scat- 
tering matrix: 



(4) 


A simple algebraic transformation leads to relations between the submatrices of 
S/ 1 and the submatrices of Pj 11 in Equation (2) 


r 1 = -p; 1 p 3 


R n = p 2 p-i 


T 1 


P - P P ~ 1 

1 r 2 '4 


T n = p-i 
4 


where 



(5) 
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The elements of S/ 1 completely describe the behavior of region II with respect 
to plane waves. It is e.g., 


R 1 




where the elements 


3 R G 


I - 




H^H 




( 6 ) 


are related to the reflection of gravity waves (G) or heat conduction waves (H) 
coming from below. 



g ^ 1 

a G 



are conversion coefficients transferring part of the incoming wave energy of one 
wave type into reflected wave energy of the other wave type. Equivalent relations 
hold for the other three submatrices R 11 , T 1 and T 11 . 

Figure 1 shows the four elements of the first column of matrix S* 1 (Equa- 
tion (4)] which give the ratios between reflected and transmitted waves, respec- 
tively, and an incoming gravity wave from below. 
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Twelve equivalent elements for the three other incoming waves exist but are 
only outlined in Figure 1. 

The elements a. and b i have the dimension /Energy/Area . ' 

Therefore R- R* and T • T* (a star indicates conjugate complex values) are 
measures of the reflected and transmitted wave energy per area. But note the 
difficulties arising from an incomplete definition of energy in Appendix E. 

3. Atmospheric Model and Numerical Techniques 

For the model atmosphere in region II we chose the Harris-Priester 
model 5 at 12 00 local time (CERA, 1965). The lower boundary of region II is 
Zj = 120 km. Below Z x an isothermal atmosphere with the temperature at the 
height Zj is added as region I. The upper boundary Z I: serves as parameter 
running from Z r to 600 km. Region III then is an isothermal atmosphere with 
the temperature in the height Z ir . For some calculations (results shown in 
Figure 7) a height Z lz = 700 km has been fixed and Z l was a parameter. 

Note that in order to make regions I and III homogeneous with respect to 
characteristic waves the coefficient of heat conductivity k in these regions must 
decrease exponentially with altitude like the pressure [see Equation (A-4)] . As 
will be seen later this is not a serious limitation for the acceptance of that 
model because below 100 km heat conduction waves are purely evanescent waves 
(their wave amplitude is attenuated by a factor larger than = 0.37 after 
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propagating 1 km in the vertical direction), and gravity wave propagation is only 
slightly influenced by heat conductivity. The region above 600 km otherwise has 
no influence on the reflection characteristics of gravity waves. 

Region II is approximated by a number of homogeneous isothermal slabs of 

thickness A Z v = 1 km. Again in each slab it is * = const. Repeated calculations 

P 

with changing AZ v 2 100 m did not show significant differences in the results 
and indicate that our model with slab thickness AZ v = 1 km is a sufficient ap- 
proximation of a realistic atmosphere. 

If one starts calculations at heights below 200 km where the attenuation rate 
of heat condition waves increases rapidly, one is confronted with a practical 
difficulty in numerical techniques. The elements of Pj 11 in Equation (1) can 
become very large in magnitude, and the elements of T 1 in Equation (5) are 
differences of two large values. If these differences are smaller than the accu- 
racy of the computer, the numerical results of the elements of T 1 become 
meaningless. In order to determine T 1 under these circumstances one has to 
repeat the calculations starting from the top of region II and calculating down- 
ward. Then the error of T 1 is small enough while now T 11 becomes meaningless 
below a certain altitude. 

4. Numerical Results 

In this section we shall give numerical results for the elements of the scat- 
tering matrix determined from Equations (D-6) and (5). 
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Figure 2 gives the magnitude of the four elements G Rj , H Rj, g tJ and H Tj 
which describe reflection, conversion, transmission and coupling into heat con- 
duction waves of incoming gravity waves from below. The frequency used is 

2 77 

co = 10” 3 sec” 1 (equivalent to a period of t = = 105 min), the horizontal wave 

number used is k x = 10“ 2 km -1 (equivalent to a horizontal wavelength of 
2 v 

k x = •£— = 628 km). 

X 

The transmission coefficient G T G decreases in this altitude range by a 
factor of 170. The mean pressure decreases in the same height interval by a 
factor of 420. Therefore the amplitude ratio of the gravity wave between 120 and 
600 km is [see Equation (A-3)] 

g T g I = °.25 

large enough to be detectable. 

The reflection coefficient G R G remains constant at altitudes above 250 km. 

Its value of 8 x 10" 2 is equivalent to a reflected energy of 6.4 x 10 -3 of the 
incoming energy. The conversion coefficient H Rj which converts part of the 
incoming gravity wave energy into reflected heat conduction wave energy, has 
the value 5 x 10 -3 and remains constant above an altitude of 200 km. The con- 
stancy of the reflection factors indicates that the region above 250 km has no 
influence on the reflection characteristics of the atmosphere below 120 km. 
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The coupling coefficient „Tj which is responsible for the transfer of energy 
from the gravity mode into heat conduction mode never exceeds 3% of the magni- 
tude of the transmission coefficient G Tj . Therefore not more than l%oof the 
gravity wave energy is coupled into heat conduction wave energy. 

In Figure 3 the transmission coefficient G T G - T w has been compared with 
the transmission coefficient T R derived from a simple ray approximation [ Equa- 
tion (D— 7) 1 . Here magnitude (Figure 3a) and phase (Figure 3b) are plotted 
versus height using the same frequency as in Figure 2 and two different values 
of k x . As seen in Figure 3 ray calculation (dashed lines) is an excellent approxi- 
mation in this range of frequency and horizontal number in heights above 200 km. 
The discrepancy between full wave and ray calculations at k x = 2.15 x 10" 3 km 
arises mainly in the height region below 200 km where the temperature gradient 
of the atmospheric model is largest. 

Note that the thickness of the whole thermosphere in terms of the vertical 
wavelength of gravity waves is only of the order of one wavelength. 


In order to test the range of validity of the ray approximation the ratio between 

T 

1 w 

the magnitudes of both transmission coefficients ~r and the difference of their 

a r 

phase values in 600 km height have been plotted versus horizontal wave number 
k x in Figure 4a again for co = 10“ 3 sec -1 . Small horizontal wave numbers 
mean steep incidence of the waves. 
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In the range k x £ 10" 3 km _1 the ray approximation breaks down. This is the 
range of horizontal wavelengths where gravity waves are totally reflected within 
the lower atmosphere (see Figure 4b). In the range k x > 10 -3 km -3 maximum 
discrepancy occurs in the phase at k x = 2.15 x 10~ 3 km -1 . But from Figure 3b 
we observe that this discrepancy already is generated within the height region 
Z < 200 km. Therefore a ray approximation is valid at this value of k x in heights 
Z £ 200 km. 

In Figure 4b the magnitudes of the four elements G R G , h R^ , G Tj and H Tj 
discussed in Figure 2 have been calculated for a region II range from 120 to 
600 km and are plotted versus horizontal wave number k x . In the range k x < 

2 x 10 _3 km" 1 the reflection coefficient | G Rj| becomes larger than 1 indicating 
the difficulty in the definition of energy (see Appendix E). Therefore the elements 
are normalized according to Equation (E-8) in this region (dashed lines) and only 
give upper limits. 

There is a maximum in the magnitude of the transmission coefficient qTq 1 near 
k x = 2 x 10" 3 km -1 showing that the atmosphere behaves like a selective trans- 
former depending on frequency and horizontal wave number. 

Otherwise the pressure amplitude of gravity waves at frequency a> = 

10“ 3 sec -1 and horizontal wave number's k x > 1.5 x 10“ 2 km -1 has dropped to 
values < 10" 1 after travelling through region II from 120 to 600 km height and 
becomes difficult to detect. 
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In the whole range of k x the coupling coefficient „Tj is negligibly small 
compared with g Tq. The conversion coefficient H R,J is also small. Gravity 
waves therefore predominate throughout the entire thermosphere, and heat con- 
duction waves are unimportant for the transport of wave energy within the 
thermosphere . 

Figures 5 and 6 calculated for the two angular frequencies co = 10“ 4 sec -1 
and co = 10“ 2 sec -1 give similar results to Figure 4. In Figure 5 we observe 
essentially the same behavior as in Figure 4 only shifted toward one smaller 
magnitude of k x . Now ray approximation is valid for horizontal wavenumbers 
k x > 10" 4 km -1 , and the maximum discrepancy in phase at k x ~ 2 x 10~ 4 km -1 
in Figure 5a is again generated below Z = 200 km. Detectable wave amplitudes 
of gravity waves are expected to have wave numbers between 2 < 10“ 4 < k x < 

2 x 10“ 3 km -1 . 

Figure 6 has been calculated for an angular frequency of co = 10“ 2 sec -1 
and a height region between 200 and 600 km because this frequency is strongly 
reflected below 200 km. Now there is no maximum in the transmission coeffi- 
cient c Tj in Figure 6b because at this frequency the wave with index G behaves 
like a gravity wave at k x > 10“ 2 km" 1 [real part of eigenvalue a G negative for 
upgoing waves (see Appendix B)] and behaves like acoustic wave at k x < 10“ 2 km -1 
(a G positive for upgoing wave). In the gravity wave range a ray approximation is 
valid. Within the acoustic range the transmission coefficient qTJ is nearly 
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constant. Coupling into heat conduction waves again is negligible within the whole 
range of k x at this frequency. 

Figure 7 (full lines) finally show reflection (hRJ 1 ), transmission (hTJ 1 ), 
conversion (qR^ 1 ) and coupling ( G T H H ) into gravity waves of downgoing heat con- 
duction waves at vertical incidence (k x = 0) and for two different frequencies 
(o)= 10~ 3 sec -1 in Figure 7a; u> = 10 -4 sec -1 in Figure 7b). The model region II 
starts now in 700 km, and the elements are calculated as functions of Zj and are 
plotted versus altitude. Here we note that some elements of the scattering 
matrix exceed the value one. This peculiar behavior has already been observed 
for ascending gravity waves at steep incidence and shows the incompleteness of 
the energy definition (see Appendix E). No attempt has been made however to 
normalize the numerical values in Figure 7 with the help of a condition equivalent 
to Equation (E-8). Therefore only the relative relations between the different 
elements can be discussed. We see from Figure 7a that in altitudes below 
500 km most of the energy of heat conduction waves of frequency o> = 10" 3 sec -1 
has been converted into gravity waves the larger part reflected, the rest trans- 
mitted into the lower regions. 

Figure 7b shows the same data but for the smaller frequency o> = 10“ 4 sec -1 . 
There below 600 km most of the energy is already converted and coupled into 
gravity waves. But now coupling predominates. Ray approximation, of course, 
becomes senseless under these conditions because of the strong coupling effect. 
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Therefore the statements in the paper of Volland (1967), which are based on ray 
approximation of vertically downward propagating heat conduction waves, must 
be modified. 


For comparison the transmission coefficient G T G n of downward propagating 
gravity waves is plotted as dashed lines in Figure 7. In Figure 7b this coefficient 
below 500 km is smaller by two orders of magnitudes than the coupling coeffi- 
cient G T H n . Therefore we expect from this result that the vertical downward 
transport of energy is mainly due to this coupling from heat conduction waves 
into gravity waves. 


Unfortunately it is not possible at the present state to derive from Figure 7 
a quantitative value for the amount of energy transfer through the thermosphere. 
We only can find from Equation (A-3) an upper limit for the ratio between the 
amplitude of a gravity wave at 120 km and the amplitude of a heat conduction wave 
at 700 km taking IqTjJ 1 | < 1 as 


c ii Pii 


| G T H n | < 1.5 x 10~ 2 


C i Pi 

which shows that these waves are no longer detectable in 120 km altitude. 


5. Concluding Remarks 

Numerical full wave calculations have been performed within the frequency 
range of gravity waves (10~ 4 ^ co < 10" 2 sec" 1 ) in thermospheric heights between 
120 and 700 km. In these heights coupling occurs between gravity waves and 
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heat conduction waves. These calculations have been compared with those for 
simple ray approximation. It has been shown that 

1. Ray treatment is an excellent approximation for obliquely upward 

propagating gravity waves at horizontal wave numbers k x £ a o. 

\ 

2. The rate of coupling from gravity wave energy into heat conduction 
energy is negligibly small for upward propagating gravity waves. 
Gravity waves therefore predominate the upward transport of wave 
energy within the thermosphere. 

3. The thermosphere reacts like a selective filter with respect to gravity 
wave. The range of optimal transmission is near k x ~ 2a). 

4. At w = NT 3 sec -1 the energy of vertically downward propagating heat 
conduction waves is mainly converted into ascending gravity waves. At 
oj = 10“ 4 sec -1 however the larger part of the energy of vertically 
downward propagating heat conduction waves is coupled into downward 
propagating gravity waves. 

These calculations suffer from one uncertainty. This is the incomplete 
definition of the energy flux for the single waves resulting in reflection and 
transmission coefficients with magnitudes larger than 1, which happens at steep 
and vertical incidence. No way has been found to overcome this difficulty. 
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APPENDIX A 


The Basic Equations 

Since we deal with gravity waves of periods 10~ 4 < a> < 10~ 2 sec" 1 at 
heights below 600 km we can neglect the coefficient of viscosity, ion drag and 
Coriolis force. We moreover consider only the simplest wave form, namely 
plane, nonducted, harmonic waves of angular frequency co propagating obliquely 
into a quiet horizontally stratified, isothermal atmosphere. The parameters of 
the atmosphere are only functions of altitude Z. The propagation plane shall be 
the (X - Z) -plane of a Cartesian coordinate system. A strict perturbation 
method then leads to a system of first order differential equations (Volland, 1968) 

% = jk Ke (A-l) 

where 



e 


l 



e 


2 
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Au 

Aw 

Ap L 
AT 

A P ^ 

cT 

p(Z) 

_ > 

T 

P(Z)> 


K 




AT 


k dAT 

e - " /|c <2 


Au 



e 


2 


A P 



~ e 3> 


= wave amplitudes of 


horizontal velocity 
vertical velocity 
< pressure 
temperature 
density 


= time independent values of 


r 


L 


acoustic phase velocity 
pressure 
temperature 
density 


coefficient of heat conductivity 

ratio of specific heats at constant pressure and constant volumn 
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k x = horizontal wave number in X-direction considered as real and constant 

00 

k = -£ wave number of acoustic waves 

co = angular frequency considered as real and constant 


C 



acoustic phase velocity 


R = gas constant 
M = molecular weight 
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_ ZI 

"a 2C 

yg 

"h 2V 

V = heat conduction velocity 

C P P 

g = gravitational acceleration 

The wave amplitudes of horizontal velocity Au and of density Ap are linearly 
related with e 2 and e 3 and do not appear explicitly in system Equation (A-l). 


The e k are normalized in such a manner that their dimension is |/Energy/Area . 
Thus the amplitude increase of the wave within an adiabatic isothermal atmos- 
phere like 



(A-2) 


(H = Tj — is the scale height) 

Mg 

is already compensated, and the squares of the normalized wave amplitudes e K e K 
are related to the energy flux of the wave. 


The ratio of the pressure amplitudes of the perturbation relative to the mean 
atmospheric pressure in two different heights Z l and Z t t is 


Pn ^Pi _ )/pi e 2^ 1 
Pi ^Ph j/Pii^n e 2 x 
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The same expressions hold for the other physical parameter. 

The height dependence of K in Equation (A-l) within an isothermal atmos- 
phere comes only through the parameter G (Z) because k = const, but 

P = p 0 e" z/H (A-4) 

Since we shall approximate the real atmosphere by a number of homogeneous 
isothermal slabs we have to choose within one slab G = const, which implies 

P «K = *0 e " Z/H (A-5) 
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APPENDIX B 


Characteristic Waves and Eigenvalues 

If the coefficient matrix K in Equation (A-l) is constant we can define as- 
cending and descending characteristic waves by 


a. 


b. 


y ~ e Tjk<,iZ 


(B-l) 


We call q A the eigenvalue of the i-th characteristic wave. Within the isothermal 
atmosphere it is 


% 



S 2 - jG ±l/(-£ 




+ 2j G(1 + B 2 S 2 ) 


(B-2) 


with 


B = 



w b 


U) 


(a> b = Brunt-Vaisala-frequency) 


The plus sign inside the square root in Equation (B-2) is associated with the pair 
of acoustic -gravity waves (an ascending one (a x ) and a descending one (b t )] . 

The minus sign is associated with heat conduction waves (a 2 , b 2 ).* 


*For reasons of clarity we shall refer in Sections 2 to 4 to a t - a G - b G 

3 2 a H “* 
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The terms "ascending” and "descending" are related to the direction of 
energy propagation. Wave energy is dissipated by heat conduction, governed by 
the parameter G in Equation (B-2), which leads to complex eigenvalues 

<Ji = «i ” J P- x ( B “ 3 ) 

From Equation (B-l) it follows then that wave energy of an upgoing wave [minus 
sign in Equation (B-l)] is dissipated if 

> 0 (B-3) 

regardless of whether a i is positive or negative. a l is negative for the 
gravity mode. This is a distinctive mark for that wave type. Within a nondis- 
sipative atmosphere G - 00 . Then acoustic -gravity waves become inattenuated 
(fiy = 0) while heat conduction waves become evanescent waves (/ 3 2 - ®). 
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APPENDIX C 


The Transformation Matrix Q 

The transformation from the physical parameters e k in Equation (A-l) into 
the parameters of the single characteristic waves a i is defined by 

e = Q a (C-l) 


where 



is the column matrix of the characteristic waves and Q can be found from the 
system of linear equations 

(K + qi E) Q a = 0 (C-2) 


where Q* is the i-th column matrix of Q and E is the unit matrix. We obtain as 
elements of Q 


F. 


Qk i ^k i 


/2qJ~ 


(C-3) 
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with 


L 


ii 



L 2 . = y- 2 A 2 - 2jAq. 

L 3i = 7 “ S 2 - A 2 - q 2 


L 4i = “2j^ ( q i + j A > 


dL 3i 


r 




1 i = 1,2 

for 

-1 i = 3,4 


q 3 ~ -c Ji 

q 4 = ~q 2 - 


F. are arbitrary normalization factors which can be taken e.g., as 


F. = 1. 

1 


But note their significance in the WKB -approximation [see Equation (D-10)] . 


Furthermore we need the reciprocal Q" 1 of the matrix Q which has the 
elements: 
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(Q 1 ) i 2 E i E 3 i + 1 E li + 2 ^ 

(Q ~ 1 )i 3 = E i CQi - j A) S i 

, 2 j G 

(Q _1 )i4 = E i~5- S i 


with 


E 


1 

0- 3i 


1 

A = 

4 S 2 A 2 
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APPENDIX D 

Exact Solution of Equation (A-l) and Ray Approximation 

Within a homogeneous slab Equation (A-l) can be transformed with help of 
Equation (C-l) into 

do 

^ = j k N o (D-l) 

with 



Within the v-th slab of thickness AZ v = Z v - Z v-1 Equation (D-l) has the solution 

« V (Z V ) = e ikvNl '^ Zl ' a v (Z J/ _ 1 ) (D-3) 

which is equivalent to Equation (B-l). 

At the boundary between two slabs total pressure, vertical velocity, tempera- 
ture and heat flux of the wave amplitude must be continuous (see Appendix E). 

This leads to (see Figure D-l) 

r e = r e (D-4) 

1 v+ 1 e v+ 1 1 v j 

Z„ + 0 Z v ~ 0 
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with 


r„(Z) = 


/ C vP, 


c y 0 0 


0 p(Z) 0 


0 0 


0 0 



(D-5) 


The general solution of Equation (A-l) for wave propagation within an atmosphere 
composed of n homogeneous slabs is then: 


%< z „> = r;> o< z o) 


(D-6) 


with 


p; = n e * ki ' N v^ z v Q -1 r;1 (Z^^) r,., (VO 0,-, 


v=l 


If the parameters of two adjacent slabs only slightly differ from each other then 


r; 1 (von.-, = . 


and the approximate ray solution of Equation (D-6) is 

i It k N d £ 

° « 0 (Z 0 ) (D-7) 

where now k(Z) q A (Z) are considered as continuous functions of altitude. The 
physical parameters of the i-th characteristic wave according to Equation (C-l) are 
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e k (Z) = Q kl (Z)a.(Z) 


(D-8) 


The approximation Equation (D-7) becomes optimal if one normalizes the elements 
of Q in such a manner that 


Q -1 


dQ 

dZ 


\ 


0 



0 


(D-9) 


contains only zeros within the diagonal. 


This condition leads to a unique determination of the factor F. in Equa- 
tion (C-3): 



(D-10) 


with 


D. 



Equation (D-8) with the normalization condition Equation (D-9) is known as WKB- 
approximation. 
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APPENDIX E 


Energy Flux 

The equation of conservation of energy neglecting viscosity can be written 
as (see e.g., Volland, 1967) 

3 

JT (E + U + X) + div S = Q (E-l) 

p 

E = — v 2 (kinetic energy) 

U = c v/ oT (internal energy) 

X = -p g(Z - Z o ) (potential energy) 

Q = (external energy source) 

S m = (E + u + X + p) V m - K (E-2) 

m 

(component of energy flux vector S in m -direction) 

Z Q is the height of an arbitrary reference level of potential energy. Its sig- 
nificance is not quite clear, and we set x = 0. The kinetic energy E is of third 
order and can be neglected. 

If there is no external energy source (Q = 0) as we assumed throughout this 
paper then the vertical component S z must be continuous at any level. This im- 
plies pressure p , temperature T, vertical velocity v z and vertical heat flux 
3 T 

- k to be continuous at any boundary. 

O Z 
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The effective energy flux of the i-th characteristic wave is given by the time 
average of S m : 



J_ y 

4 (7 - 1) 


2S 


e e* 
e 2 e 2 


(e 


e * + 
2 e l 


e i) 


(E-3) 


(The star indicates complex conjugate values.) The vertical time averaged heat 
flux is canceled out because 


K 


d AT 
dZ 



e 


4 


in Equation (A-l) and 


e „ = 0 


Unfortunately the direction 6 of energy flux derived from Equation (E-3) 


tg 6 

* e 1 S, 


S [(7 - 2 A 2 - 2 A/3j ) 2 + 4 A 2 a 2 ] 
a. (y 2 - 4 A 2 S 2 ) 


(E-4) 


does not coincide with ray direction 6 1 derived by 


t R0. 


d a i 

~Ts 


= Real J 




1 - 


GB 2 


q A 2 + A 2 - ~ + S 2 + j G 


(E-5) 
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as long as k ^ 0. Only for G -> oo (/< - 0) do Doth directions become equal for 
gravity waves: 

tg^ e = tg0 r = - (1 - B2) (* = 0) (E-6) 

We observe the same discrepancy as in electromagnetic wave propagation where 
the direction of the Poynting vector and the ray direction differ from each other 
within a dissipative plasma. (See e.g., Hines, 1951.) Equation (E-3) therefore 
is no definite measure of energy flux. Whether or not the neglect of potential 
energy is the reason for this discrepancy is not clear to us. 

In view of these difficulties we do not try to use Equation (E-3) for the de- 
termination of energy flux but simply consider the values a i a* calculated from 
solution Equation (D-6) as measure of energy flux of the different characteristic 
waves. 

If the amplitudes a. are clearly related to wave energy flux then the magni- 
tudes of the elements of the scattering matrix [see Equation (4)] should be always 
smaller than 1. In the numerical results of Section 4 we shall find, especially 
at steep incidence, magnitudes of some reflection and transmission coefficients 
larger than 1. This indicates that our definition of energy flux is incomplete 
resulting from the difficulties outlined above. 
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From general energy considerations however we can give an upper limit 
for the magnitude of the elements of the scattering matrix. Imagine that only a 
gravity wave penetrates into region II from below (see Figure 1). Then Equa- 


tion (3) leads to 



(E-7) 


Since all waves are decoupled from each other in regions I and III and since 
energy leaving region II can not exceed energy penetrating region II then 

(g R g ) 2 + (hRq 1 ) 2 + (qTg 1 ) 2 + (hTJ ) 2 ^ 1 (E-8) 

where the equal sign in Equation (E-8) applies only in the case of a nondissipative 
atmosphere (k = 0). 

From Equation (E-8) and equivalent equations holding for the other 12 ele- 
ments we can normalize the R and T. This has been done partly in Figures 4 
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FIGURE CAPTIONS 

Figure 1. Elements of the scattering matrix. Reflection coefficients (R) and 
transmission coefficients (T) of upward (a) and downward (b) propagating 
gravity waves (G) and heat conduction waves (H) within the thermosphere. 

Figure 2. Magnitude of coefficients of reflection (qRJ), transmission ( G Tj ), 
conversion ( H Rj ) and coupling ( H Tj) into heat conduction waves of upward 
propagating gravity waves versus Z ir . 

Figure 3. Magnitude (Figure 3a) and phase (Figure 3b) of transmission coeffi- 
cients of upward propagating gravity waves calculated by full wave theory 
(full lines) and by ray approximation (dashed lines). 

Figure 4a. Ratio between transmission coefficients of upward propagating 
gravity waves calculated by full wave theory (W) and ray theory (R) versus 
horizontal wave number k x . Range of thermospheric model: 120 to 600 km; 
angular frequency: a > = 10" 3 sec 1 . 

Figure 4b. Magnitude of coefficients of reflection, transmission, conversion 

and coupling into heat conduction waves of upward propagating gravity waves 
versus horizontal wave number k x . Range of thermospheric model: 120 to 
600 km. Angular frequency: oj = 10" 3 sec' 1 . Dashed lines give upper limits 


(for details see text). 
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Figure 5. Same as text of Figure 4, but angular frequency: = 10' 4 sec' 1 . 

Figure 6. Same as text of Figure 4, but angular frequency: co = 10' 2 sec' 1 ; 
range of thermospheric model 200 - 600 km. 

Figure 7. Full lines: Magnitude of coefficients of reflection, transmission, 
conversion and coupling into gravity waves of vertically downward propa- 
gating heat conduction waves versus Z I . 

Dashed lines: Magnitude of transmission coefficients of vertically 
downward propagating gravity waves . 

Figure 7a. co = 10' 3 sec* 1 . 

Figure 7b. co = 10' 4 sec" 1 . 

Figure D-l. The homogeneous slab-model. 




Figure 1. Elements of the scattering matrix. Reflection coefficients (R) and transmission co- 
efficients (T) of upward (a) and downward (b) propagating gravity waves (G) and heat conduction 
waves (H) within the thermosphere. 


37 












HORIZONTAL WAVE NUMBER k^km' 1 ) 



10 4 10 3 10 2 


ll/A 



HORIZONTAL WAVE LENGTH (km) 


Figure 4a. Ratio between transmission coefficients of upward propagating gravity waves calcu- 
lated by full wave theory (W) and ray theory (R) versus horizontal wave number k x . Range of 
thermospheric model: 120 to 600 km; angular frequency: co ~ 10" 3 sec' 1 . 
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Figure 4b. Magnitude of coefficients of reflection, transmission, conversion and coupling into 
heat conduction waves of upward propagating gravity waves versus horizontal wave number k x . 
Range of thermospheric model: 120 to 600 km. Angular frequency: co - 10" 3 sec"\ Dashed 
lines give upper limits (for details see text). 
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Figure 5b. Same as text of Figure 4, but angular frequency: co - 10- 4 sec" 1 . 
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Figure 7a. Full lines: Magnitude of coefficients of reflection, transmission, conversion and 

coupling into gravity waves of vertically downward propagating heat conduction waves versus Zj. 

Dashed lines: Magnitude of transmission coefficients of vertically downward propagating 

gravity waves. Angular frequency: co - 10‘ 3 sec'*. 
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Figure 7b. Full lines: Magnitude of coefficients of reflection, transmission, conversion and 
coupling into gravity waves of vertically downward propagating heat conduction waves versus Zj. 

Dashed lines: Magnitude of transmission coefficients of vertically downward propagating 

gravity waves. Angular frequency: co - 10’^ sec"^. 
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